Load required libraries

library(prettydoc)
library(tidyverse)
library(gsw)
library(readr)

Now we need to import our data

hydrostation_bottle <- read_delim("hydrostation_bottle.txt", 
    delim = "\t", escape_double = FALSE, 
    col_names = FALSE, trim_ws = TRUE, skip = 31)

hydrostation_bottle_names <- read_csv("hydrostation_bottle.txt", 
    skip = 30)

colnames(hydrostation_bottle) = colnames(hydrostation_bottle_names)
#view(hydrostation_bottle)

Hydrostation S Discrete Bottle Data for years 1955 through December 2020.

Variable Names and Units

yyyymmdd = Year Month Day
decy = Decimal Year
time = Time (hhmm)
latN = Latitude (Deg N)
lonW = Longitude (Deg W)
Depth = Depth (m)
-Temp = Temperature ITS-90 (C)
Pres = CTD Pressure (dbar)
CTD_S = CTD Salinity (PSS-78)
Sal1 = Salinity-1 (PSS-78)
Sig-th = Sigma-Theta (kg/m^3)
O2(1) = Oxygen-1 (umol/kg)
OxFixT = Oxygen Fix Temp (C)
Anom1 = Oxy Anomaly-1 (umol/kg)
Quality flags
-999 = No data
0 = Less than detection limit

# lets first plot the data
hydrostation_bottle %>% 
  filter(`Sig-th`!= -999) %>% # filter out -999 no data flag
  ggplot()+geom_point(aes(x=decy, y = `Sig-th`)) # hard to interpret even w no -999s

hydrostation_bottle %>% 
  filter(`Sig-th`!= -999 & Depth <20) %>% # filter out -999 no data flag and by upper 20m
  ggplot()+geom_line(aes(x=decy, y = `Sig-th`)) # line shows seasonality

# clear seasonal signal for sigma-theta, lets see how this compares to temp
hydrostation_bottle %>% 
  filter(`Sig-th`!= -999 & Depth <20) %>% # filter out -999 no data flag and by upper 20m
  ggplot()+geom_point(aes(x=Temp, y = `Sig-th`))

TEOS-10 Toolbox in Package Seacarb

?gsw #launches documentation for gibbs seawater toolbox

## starting httpd help server ... done

? gsw_sigma0 # lets check this function 
# it says we need absolute salinity and conservative temperature

#first we need absolute salinity

?gsw_SA_from_SP
#practical salinity
#sea pressure(dbar)
#longitude
#latitude

#plot our pressure data - its missing before 1980s
hydrostation_bottle %>% 
  ggplot()+geom_point(aes(x=decy,y=Pres))

#we have depth for the time series
hydrostation_bottle %>% 
  ggplot()+
  geom_point(aes(x=decy, y=Depth))

#adds a pressure column from depth and latN columns from/to hydrostation bottle
hydrostation_bottle = 
  hydrostation_bottle %>% 
  mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,))

hydrostation_bottle %>% 
  ggplot()+geom_point(aes(x=Pres,y=Pres_gsw))

## Warning: Removed 3 rows containing missing values (`geom_point()`).

# we see strong 1:1 agreement between measured pressure and calculated pressure

hydrostation_bottle %>% 
  ggplot()+geom_point(aes(x=decy,y=latN))

hydrostation_bottle = 
  hydrostation_bottle %>% 
  mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>% 
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN))
#plot it
hydrostation_bottle %>% 
  filter(Sal1!=-999) %>% 
  ggplot()+
  geom_point(aes(x=Sal1,y=S_abs_gsw))

# now we need to calculate conservative temperature
# we need absolute salinity, insitu temp (ITS-90), and sea pressure
hydroS =
  hydrostation_bottle %>% 
  filter(Sal1!=-999) %>%
  mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>% 
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>% 
  mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw))

hydroS %>% 
  filter(Temp!= -999) %>% 
  ggplot()+
  geom_point(aes(x=Temp, y =T_cons_gsw))

#add line to calculate conservative temperature
hydroS =
  hydrostation_bottle %>% 
  filter(Sal1!=-999) %>%
  mutate(Pres_gsw = gsw_p_from_z(Depth*-1,latN,)) %>% 
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>% 
  mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>% 
  mutate(Sig_th_gsw = gsw_sigma0(S_abs_gsw,T_cons_gsw))

hydroS %>% 
  filter(`Sig-th`!= -999) %>% 
  ggplot()+
  geom_point(aes(x=`Sig-th`, y=Sig_th_gsw))

hydroS %>% 
  filter(Sig_th_gsw<0)

## # A tibble: 544 × 19
##         Id yyyymmd  decy  time  latN  lonW Depth  Pres  Temp CTD_S  Sal1 Sig-t…¹
##      <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
##  1  6.00e9  1.96e7 1955.  1703  32.2  64.4   495  -999  -999  -999  36.3    -999
##  2  6.00e9  1.96e7 1955.  1601  32.1  64.3   144  -999  -999  -999  36.6    -999
##  3  6.00e9  1.96e7 1956.  1704  32.1  64.3   873  -999  -999  -999  35.4    -999
##  4  6.00e9  1.96e7 1956.  1704  32.1  64.3    50  -999  -999  -999  36.2    -999
##  5  6.00e9  1.96e7 1956.  1704  32.1  64.3   244  -999  -999  -999  36.5    -999
##  6  6.00e9  1.96e7 1956.  1704  32.1  64.3   685  -999  -999  -999  36.0    -999
##  7  6.00e9  1.96e7 1956.  1705  32.2  64.3    50  -999  -999  -999  36.4    -999
##  8  6.00e9  1.96e7 1957.  1703  32.1  64.3  2024  -999  -999  -999  35.0    -999
##  9  6.00e9  1.96e7 1958.  1504  32.0  64.4   961  -999  -999  -999  35.2    -999
## 10  6.01e9  1.96e7 1958.  1705  32.2  64.2   553  -999  -999  -999  36.1    -999
## # … with 534 more rows, 7 more variables: `O2(1)` <dbl>, OxFix <dbl>,
## #   Anom1 <dbl>, Pres_gsw <dbl>, S_abs_gsw <dbl>, T_cons_gsw <dbl>,
## #   Sig_th_gsw <dbl>, and abbreviated variable name ¹`Sig-th`

  #view()


# how to replace data in a line without making this many dataframes (i.e. select and/or filters)
hydroS_correctedS_a =
  hydroS %>% 
  filter(Sig_th_gsw<0) %>%
  mutate(S_abs_gsw=gsw_SA_from_SP(CTD_S,Pres_gsw,360-lonW,latN)) %>% 
  mutate(T_cons_gsw = gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>% 
  mutate(Sig_th_gsw = gsw_sigma0(S_abs_gsw,T_cons_gsw)) 


hydroS_correctedS_b =
  hydroS %>% 
  filter(Sig_th_gsw>0)

hydroS_corrected = rbind(hydroS_correctedS_a, hydroS_correctedS_b)

hydroS_corrected %>% 
  filter(`Sig-th`!= -999) %>% 
  ggplot()+
  geom_point(aes(x=`Sig-th`, y=Sig_th_gsw))

## Next day
hydroS_corrected %>% 
  ggplot()+
  geom_point(aes(x=Sig_th_gsw, y = Depth))+
  scale_y_reverse()+
  scale_x_continuous(position ="top")+
  xlab(expression(paste(sigma[theta]," (kg m"^"-3",")")))+
  ylab("Depth(m)")

## Warning: Removed 543 rows containing missing values (`geom_point()`).

Has surface sigma theta decreased over time?

hydroS_shallow = hydroS_corrected %>% 
  filter(Depth<30)

sig_theta_tm = lm(Sig_th_gsw~decy,data =  hydroS_shallow) # y (Sig) is a function of x(decy), lm(y~x,data = data )
summary(sig_theta_tm)

## 
## Call:
## lm(formula = Sig_th_gsw ~ decy, data = hydroS_shallow)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5789 -0.8683  0.1070  0.8819  1.5805 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.4047230  1.6704055  19.998  < 2e-16 ***
## decy        -0.0042378  0.0008387  -5.053 4.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9401 on 3410 degrees of freedom
##   (56 observations deleted due to missingness)
## Multiple R-squared:  0.007431,   Adjusted R-squared:  0.00714 
## F-statistic: 25.53 on 1 and 3410 DF,  p-value: 4.587e-07

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

plot = hydroS_shallow %>% 
  ggplot()+
  geom_point(aes(x=decy, y=Sig_th_gsw))+
  geom_line(aes(x=decy, y=Sig_th_gsw))+
  geom_smooth(aes(x=decy, y=Sig_th_gsw),method = "lm")+
  theme_classic()
ggplotly(plot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 56 rows containing non-finite values (`stat_smooth()`).

#y=mx + b

#y = Sig_th_gsw
#x - decy
#coeffiecients: intercept = b, decy = m
#Sig_th_gsw = -0004*decy + 33.4

#Lab Assignment 1. Pick a question (include hypothesis)

How does oxygen look like at depth and over time? H1 : Oxygen reaches a max in the upper 1000m before decreasing with further depth

Produce a plot and statistal summary using lm()

o2_h = hydroS_corrected %>% 
  filter(`O2(1)`!= -999)
o2_h = 
  ggplot()+
  geom_point(aes(x=`O2(1)`,y = Depth))

Describe your results, the summary, and answer the question
Compile into a completed lab report using R markdown

Potential Questions :

How do temperature, salinity, and sigma-theta co-vary? Is there a relationship betwen sigma-theta and oxygen? Is there a relationship of any of the parameters with depth? With time? Within a depth range over time? Are there seasonal differences in any of the parameters?

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

TEOS Hydrostation S Data Lab

Jose Martinez

2023-03-07